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1.0  INTEQDUCTION 


Ionospheric  electron  density  distributions  show  great  longitudinal  and  latitudinal 
variations  in  space,  as  well  as  diurnal,  seasonal  and  solar  cycle  variations  in  time, 
even  for  the  quiet  ionosphere.  Generally,  ionospheric  variations  appear  in  a  wave¬ 
like  form  with  the  wavelengths  ranging  from  meters  to  several  thousand 
kilometers.  The  electron  density  of  the  ionosphere  is  a  function  of  the  location 
coordinates  X  and  Y,  the  altitude  z,  and  the  time  t: 

N  =  N(X,Y,z,t)  (1) 

The  task  of  ionospheric  mapping  is  to  find  the  function  (1)  as  the  best  fit  to  the 
observed  data  from  the  sounding  stations.  Since  the  Digisonde  network/ 1/ 
produces  electron  density  profiles  in  real  time,  the  potential  exists  to  generate 
electron  density  maps  in  real  time  for  areas' with  adequately  located  Digisondes. 
This  paper  describes  the  new  three-dimensional  mapping  technique  and  shows 
preliminary  results.  After  introducing  the  physical  and  mathematical  foundation, 
the  paper  applies  the  technique  to  the  Northeastern  American  coastal  region  using 
data  of  December  1990. 

The  Digisondes /2/  specify  the  electron  density  profile  by  a  set  of  parameters/ 3,4/  for 
each  ionospheric  layer; 


^Tn>  Ao/  Ai,  A2,  A3,  A4  (2) 

and  three  of  such  sets  determine  the  profiles  for  the  E,  the  FI  and  the  F2  layers 
through  the  relation: 

4 

z  =  Zm  +  Vi  5^  Ai  T(g)  (3) 

i=0 


with: 

ln(f/fnO 

S-ln(fs/U 


(4) 


1 


» 

where  z  denotes  the  altitude  and  f  the  plasma  frequency.  T.  is  the  i-th  order  shifted 

Chebyshev  polynomial/5/.  The  subscripts  s  and  m  denote  the  values  at  the  starting 
height  and  the  peak  of  the  layer,  respectively.  For  the  F  region  mapping,  we  find  a 
new  function  (3),  whenever  the  FI  layer  exists,  which  describes  the  profile  over  the 
entire  F  region/6/.  This  procedure  avoids  the  discontinuity  problem,  caused  by  the 
appearance  of  FI  at  some  but  not  all  stations  in  the  mapping  region,  and  causes  only 
minimal  distortions. 

The  general  problem  is  to  determine  the  so-called  objective  function  for  the 
parameter  p: 

p  =  P(X,Y,t)  (5) 

based  on  the  observed  data 
Pmn  =  P(XnvYnvfn) 

(m  =  1,  2, M;  n  =  -N+1, ...,  -2,  -1,  0) 

where  M  specifies  the  number  of  stations  involved,  to  denotes  the  present  universal 
time  and  N  is  the  number  of  samples  going  backward  in  time.  For  a  given  station 
the  data  samples  form  a  time  series  with  the  interval  x,  and  the  available  data  base 
consists  of  M  sets  of  time  series.  In  our  application  x  =  30  min,  N  =  64  and  Ln+i  = 
t-63  =  to  -  31 .5h.  By  nowcasting  we  mean  calculation  of  N(X,  Y,  z,  to)  for  the  current 
time  to,  as  illustrated  later  in  this  paper.  The  technique  also  lends  itself  to  short¬ 
term  forecasting,  i.e.  for  calculation  of  a  3D  map  for  t  >  to.  We  do  not  discuss  this 
application  here. 

The  above  mathematical  presentation  of  the  problem  assumes  that:  (1)  at  any 
arbitrary  location  (X,Y)  in  the  area  and  at  any  time  t,  the  electron  density  profile  can 
be  expressed  in  the  same  way  as  at  the  observation  sites;  and  (2)  adequate  data 
density  (both  spatially  and  temporally)  exists  to  ensure  proper  sampling  of  the  main 
structures  and  variations  in  the  ionosphere.  The  first  assumption  simplifies  the 
problem  from  finding  the  electron  density  profile  to  finding  the  function  p, 
equation  (5),  thus  eliminating  the  variables  z.  The  second  assumption  limits  the 


process  to  spatial  variations  with  wavelengths  larger  than  2L,  and  temporal 
variations  larger  than  2x,  where  t  is  the  sampling  interval  and  L  the  average  station 
separation. 

2.0  PHYSICAL  AND  MATHEMATICAL  FOUNDATION 

Investigations  of  the  ionospheric  formation,  motion  and  dependence  on  solar- 
terrestrial  parameters  lead  to  physical  models  of  the  ionosphere.  Because  of  the 
complicated  physical  processes  involved  the  theoretical  models  are  generally  less 
accurate  than  models  derived  from  measured  data.  Our  data  drive  technique 
derives  from  the  concept  that  changes  of  the  ionosphere  at  any  one  location  relate  to 
those  at  another  location  (at  least  within  a  limited  area),  and  that  ionospheric 
plasma  changes  occur  in  a  wave-like  form.  For  example,  the  diurnal  variation  at  a 
western  location  lags  in  time  with  respect  to  a  more  eastern  location.  For  a 
morphological  description,  the  time  history  of  the  sampled  data  contains  sufficient 
information  to  analyze  the  behavior  of  the  ionospheric  changes.  The  technique 
decomposes  the  ionospheric  variations  in  the  mapping  region  into  a  series  of 
simple  plane  waves  which  describe  the  changing  ionosphere  at  any  location. 
Fourier  analysis  provides  a  tool  for  the  wave  decomposition  and  we  assume  the 
existence  of  a  Fourier  expansion  of  function  (5): 

+00 

P(X,Y,t)  =  JJJ  P(o),kx,ky)  dto  dkx  dky  (7) 

•CO 

Equation  (7)  decomposes  the  time-varying  two-dimensional  physical  field  P(X,Y,t) 
into  a  series  of  plane  waves.  The  relation  between  the  spatial  wave  numbers  kx  and 
ky  and  the  frequency  to,  i.e.  the  dispersion  relation,  derives  from  the  dynamic 
equations  of  the  field  and  the  properties  of  the  medium.  We  do  not  intend  to  study 
the  physical  aspects  of  the  field,  but  simply  assume  that  kx  and  ky  are  complex 
functiorts  of  to  only,  i.e.: 

kx  =  kx(to)  +  iSx(to)  ,  ky  =  ky(to)  +  iSy(to)  (8) 
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This  means  that  the  parameters  of  a  wave  propagating  in  a  given  direction  do  not 
depend  on  location.  The  objective  function  becomes: 


4*00 


P(X,Y,t)=  J  P(o)) 


gSx(a))X+Sy(a))Y+i[cot-kx(a»X-ky(to)Y] 


Rewriting  Pfco)  as: 

Pto)  =  A((o)  <9) 

we  have: 


4-O0 


P(X,Y,t)=  J  A(e)) 

•oo 


gSx((o)X+Sy{a))Y+i[a)t-kx(a))X-ky((o)Y-«»o(“)) 


(10) 


where  A(o)),  Sx(to),  Sy(o)),  kx(o)),  ky(co)  and  are  all  real  functtons. 

For  any  time  tn^to,  we  calculate  the  objective  function  P(X,Y,t)  with  the  help  of  the 
discrete  Fourier  transform  (DFT).  With  catn  =  27cvn/N,  P(X,Y,tn)  becomes  /7/; 

N/2-1 

P(X,Y,tn)=  X  (11) 

V  =  -N/2 

Equation  (11)  decomposes  the  time-varying  two-dimensional  physical  field  into  N 
damped  plane  waves  with  different  frequencies  and  amplitudes.  The  only 
important  limitation  possibly  restricting  application  of  (11)  resides  in  equation  (8) 
which  assumes  that  the  parameters  of  any  wave  in  a  given  direction  do  not  change 
along  the  wave  trajectory.  In  this  way,  the  mapping  problem  reduces  to  finding  the 
functions  A(v),  Sx(v),  Sy(v),  kx(v),  ky(v)  and  Oofv).  The  next  section  will  show  how  to 
resolve  the  difficulty  of  phase  ambiguity  and  how  to  determine  these  functions. 
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3.0  DETERMINATION  OF  THE  FOURIER  KERNEL 


At  any  station  X=Xm  and  Y=Ym: 

N/2-1 

P(XnvYm/tn)=  ^  y^^yjgSxXm+SyYnTiIkxX}n+kyYu,+<t>ol+i2jtnv/N 
V  =  -N/2 


(12) 


We  expand  the  measured  data  series  (6)  also  into  a  Fourier  series: 

N/2-1 

P(XnvYnvtn)=  ^  Pm(v)e'2''"''/^  (13) 

V  =  -N/2 

Expressing  Pin(v)  in  terms  of  amplitude  and  phase 

Pm(v)  =  Pom(v)  ,  (14) 

(-It  <  <t>m(v)  <  7l)  , 

and  comparing  (12)  to  (13),  we  obtain 

Pom(v)  =  A(v)  (15) 

and 

^m(v)  =  hxXm+kyYm+<l>o+2km7l  ,  (16) 

where  kni(v)  is  an  integer,  generally  0  or  ±1. 

For  three  stations  (M  =  3)  in  the  mapping  area,  equations  (15)  and  (16)  determine 
values  of  A,  Sx,  Sy,  kx,  ky  and  <I)o  for  each  v.  For  more  than  three  stations,  a  best-fit 
technique  determines  the  values. 

The  error  functions 


5 


M 


ei  =  I  (InPom(v)  -  InA(v)  -  S^Xm  -  SyYmP 

m=l 


(17) 


M 

£2  "  ^  [0m(v)  -  2km7t  -  kxXm  -  kyYm  -  (j)o]^ 

m=l 


(18) 


will  be  minimized  if  the  A(v),  Sx,  Sy,  kx,  ky  and  0o  satisfy  the  following  two  sets  of 
relations: 


M  MM 

lnA(v)M  +  S^  I  X^  +  Sy  I  Yj„=  I  InPo^(v) 

m=l  m=l  m=l 


M  M  M  M 

InA(v)  S  +  I  X2,  +  Sy  X  =  X  lnPo„(v)X„,  (19) 

m=l  m=l  m*l  m=l 


MM  MM 

InA(v)  I  Y„  +  S,  I  X„Y„  +  Sy  I  Y?,.  I  lnPom(v)Y„ 

m=l  ni=l  m=l  m=l 

and 

M  MM 

0oM  +  k^  I  X^  +  l^y  2:  Y^=  I  [0^(v)-2k^7l] 
m=l  m=l  m=l 

M  M  M  M 

Oo  I  X„  +  k,  I  X^  +  ky  I  X„Y„=  I  |$„(v)-2k„*lx„  (20) 

m=l  m=1  m=l  m=1 

MM  MM 

♦0  X  S  V„*ky  I  Y2,=  X  14.„(v)-21s„7tlY„ 

m=l  m=l  m=l  m=l 


These  are  first-order  linear  equations.  From  the  first  set  of  equations  (19),  we  find 
A(v),  Sx  and  Sy,  and  from  the  second  set  of  equations  (20)  we  find  0O/  and  ^y- 
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Consequently,  the  objective  function  is  determined  and  the  parameter 
P(X,Y,tn)  is  calculated  with  equation  (11)  for  any  location,  assuming  we  can  estimate 
the  correct  value  for  the  integer  in  Eq.  (20).  For  the  low  frequencies,  it  is  expected 
that  the  values  of  <I)m  for  different  stations  are  close  to  each  other,  so  there  is  no 
phase  ambiguity.  For  higher  frequencies  different  sets  of  values  (0n;±1)  are  tested 
as  illustrated  in  the  following  example.  A  trial  technique,  with  the  least  error 
criteria,  generally  finds  the  right  values  for  k^.  Wrong  decisions  on  the  k^  values 
can  occur  near  the  Nyquist  frequency  l/(2x)=lh-l.  Fortunately,  these  higher  spectral 
components  with  about  1  hour  periods  in  our  analysis  correspond  to  gravity  waves 
that  contain  little  energy  and  do  not  introduce  serious  errors  in  the  resulting  maps. 
The  following  example  illustrates  the  procedure.  The  observed  phases  Om  (see  Eq. 


14)  were  <l>i  = 

12.3°,  <1)2 

=  89.8°,  <1)3 

—  “ 

113.5°,  <1)4  = 

=  -170.1 

°,<1>5  = 

■22.4. 

Case  # 

ki 

k2 

k3 

k4 

ks 

<1), 

<1)2 

(1)3 

<1)4 

<1)5 

1 

0 

0 

0 

0 

0 

12.3 

89.8 

-113.5 

-170.1 

-22.4 

2 

0 

0 

0 

-1 

0 

12.3 

89.8 

-113.5 

189.9 

-22.4 

3 

0 

0 

-1 

-1 

0 

12.3 

89.8 

246.5 

189,9 

-22.4 

4 

0 

0 

-1 

-1 

-1 

12.3 

89.8 

246.5 

189.9 

237.6 

5 

-1 

0 

-1 

-1 

-1 

372.3 

89.8 

246.5 

189.9 

237.6 

A  least-squares  fitting  is  done  for  each  case  and  the  one  with  the  best  fit  is  adopted. 

4  0  APPLICATION  OVER  THE  NORTHEASTERN  AMERICAN  COAST 

We  applied  the  new  technique  to  mapping  the  electron  density  structure  over  the 
northeast  American  coastal  region  within  the  area  30°  to  60°N  and  80°  to  50°W. 
Five  Digisonde  stations  (Figure  1)  supplied  the  data  for  the  regional  mapping  during 
December  1990; 

(53.32°N,  60.33°W) 

(47.28°N,  53.97°W) 

(42.60°N,  72.50°W) 

(37.93°N,  75.47°W) 

(32,33°N,  64.67°W) 


Goose  Bay,  Labrador 
Argentia,  NF 
Millstone  Hill,  MA 
Wallops  Island,  '^A 
Bermuda 
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Figure  1.  Five  Digisonde  Stations  at  the  northeast  American  coastal  region  (G 
Goose  Bay,  A-Argentia,  M-Millstone  Hill,  W-Wallops  Island  and  B 
Bermuda)  make  routine  half-hourly  ionograms. 
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Figure  2  shows  the  time  variations  of  one  of  the  profile  parameters,  foF2,  for  the 
five  stations.  Small  dots  mark  the  data  samples  used  to  calculate  the  spectrum  for 
1300  UT  on  day  349,  1990.  For  each  station  a  32  hour  time  series  with  30  minute 
samples  (N=64,  t=30min)  was  Fourier  transformed  using  an  FFT  algorithm, 
resulting  in  the  spectra  of  Figure  3.  The  64  line  amplitude  spectrum  on  the  left  has 
even  symmetry,  and  the  64  line  phase  spectrum  on  the  right  has  odd  symmetry  as 
expected  for  a  real  function.  In  the  plot  a  value  of  27t  had  been  added  to  negative 
phases.  The  spectral  lines  represent  the  values  of  Pom(v)  and  respectively,  in 

equation  (14).  Figure  3  illustrates  the  typical  behavior  of  all  observed  amplitude 
spectra;  a  rapid  amplitude  decrease  with  frequency.  This  means  that  the  lower 
frequencies  dominate  the  three-dimensional  profile  reconstruction. 

Figure  4  illustrates  the  wave  superposition  technique  displaying  the  time 
development  of  the  foF2  parameter  in  the  30°  by  30°  mapping  region  during  sunrise, 
from  1100  UT  to  1430  UT.  The  white  dots  indicate  the  locations  of  the  five  stations. 
The  gray  shade  scale  on  the  right  side  reveals  that  the  plasma  frequencies  decrease 
from  above  5  MHz  in  the  South-East,  to  below  4  MHz  in  the  North-West  at  1100  UT, 
and  from  above  13  MHz  to  below  11  MHz  at  1430  UT.  As  a  consistency  check  the 
measured  foF2  values  at  each  of  the  five  stations  are  shown  on  top  of  the  figure 
together  with  the  "mapped"  value  (in  parenthesis).  For  this  example  of  a  sunrise 
period,  the  differences  between  measured  and  mapped  values  are  acceptably  small. 
Figure  4  (on  top)  also  indicates  the  available  data  samples  for  each  station.  In  the 
current  analysis,  data  gaps  were  filled  by  linear  interpolations. 

Accurate  HF  ray  tracing  between  two  separate  points  requires  the  knowledge  of  the 
electron  distribution  in  a  vertical  plane  through  the  two  points.  This  requires  the 
calculation  of  vertical  profiles  along  the  great  circle  path  between  the  two  points. 
From  the  objective  functions  of  all  parameters  specified  in  (2)  the  program  calculates 
the  two  dimensional  electron  distributions.  Figure  5  shows  the  sequence  of  east- 
west  profile  cross  sections  through  the  points  (280°,  45°)  and  (300°,  45°)  during 
sunrise  on  day  90-349  in  half-hour  increments.  Starting  at  1100  UT  tilts  develop  at 
305°E  longitude  which  by  1130  UT  extend  over  the  entire  2500km  region.  At  1330 
UT  the  F  layer  is  generally  horizontally  satisfied,  except  of  a  small  ondulation. 
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Figure  2.  Time-variation  of  parameter  foF2  for  5  stations. 

To  calculate  the  contour  map  for  foF2  at  1300  UT  on  day  349,  the 
previous  32  hours  of  data  for  each  station  are  Fourier  transformed.  The 
data  samples  used  for  constructing  the  1300  UT  map  are  marked  by  small 
dots. 
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Figure  3  Spectra  for  parameter  foF2  for  5  stations  day  349, 1990, 13  UT. 

The  amplitude  spectra  determine  which  "wave  frequencies"  contribute 
to  the  parameter  map.  The  amplitude  and  phase  distribution  for  each 
spectrum  component  over  the  five  stations  specifies  the  attenuation  and 
phase  constants  and  the  wave  directions.  Linear  superposition  of  the  set 
of  plane  waves  generates  the  parameter  contour  map  for  the  30°  by  30° 
area.  Once  the  contour  maps  for  all  parameters  are  specified,  the  electron 
density  profile  at  any  point  in  the  area  is  given  by  equation  (3). 
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Figure  4.  Sunrise  sequence  of  foF2  maps. 

Systematic  change  of  the  foF2  contours  as  a  function  of  time;  half-hourl) 
increments. 


The  nowcasting  software  ran  on  a  386  PC  computer  in  the  laboratory.  In  an 
operational  environment  all  new  data  arriving  from  the  stations  can  feed  real-time 
calculations  for  continuous  generation  of  updated  maps.  This  method  of 
nowcasting  can  support  OTH  radar  operation  as  well  as  HF  communications  and 
direction  finding  applications.  The  updated  electron  density  maps  will  "follow"  the 
changes  of  the  ionosphere.  Furthermore,  using  the  measured  time  series  of  the 
parameters  from  the  monitoring  stations,  this  mapping  technique  can  be  expanded 
to  forecast  near  future  maps  of  electron  density  distribution. 

5.0  EVALUATION  OF  THE  NOWCASTING  TECHNIQUE 

Evaluation  of  the  nowcasting  technique  requires  additional  sounding  stations 
located  within  the  mapping  area,  sufficiently  separated  from  the  given  monitoring 
stations.  Data  from  one  ionosonde  in  Maine  showed  excellent  agreement  with 
mapped  values,  but  the  distance  to  the  Millstone  Hill,  MA  sounder,  which  is  one  of 
the  network  sounders,  measures  only  600  km. 

To  generate  a  broad  testing  database,  we  used  the  IRl  model/8/  to  provide  time 
sequences  of  profiles  for  the  locations  of  the  five  Digisondes  and  for  16  other 
locations  in  the  mapping  area.  The  nowcasting  algorithm  used  the  profiles  from  the 
five  sounder  locations  to  calculate  the  profiles  for  each  of  these  other  locations. 
Averaging  the  differences  between  the  database  and  mapped  values  over  a  24  hour 
period  provides  a  good  measure  of  the  performance  of  the  nowcasting  technique. 

Figure  6  shows  the  rms  errors  in  foF2  and  hmF2  for  the  months  of  January  and  July, 
both  for  sunspot  values  of  R=10  and  R=150.  The  errors  remain  small  within  the 
mapping  region,  and  even  in  the  boundary  regions  the  errors  do  not  grow 
excessively. 
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349/1990 


Figure  5.  Sunrise  sequence  of  profile  cross  sections. 


Profile  cross  sections  in  any  direction  can  be  calculated  through  the 
mapped  area.  The  time  sequence  of  the  east-west  cross  section  during 
sunrise  on  day  349  in  half-hour  increments  is  shown. 


Figure  5.  Sunrise  sequence  of  profile  cross  sections  (Concluded). 

Profile  cross  sections  in  any  direction  can  be  calculated  through  the 
mapped  area.  The  time  sequence  of  the  east-west  cross  section  during 
sunrise  on  day  349  in  half-hour  increments  is  shown. 
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Figure  6. 
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Evaluation  of  the  mapping  technique  using  the  IRI  database.  RMS  errors 
for  foF2  and  hmF2  at  the  5  sounder  sites  and  16  other  locations  for 
January  and  July  under  low  (R=10)  and  high  (R=150)  solar  activity, 
respectively. 
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6.0  SUMMARY 


The  new  data  drive  mapping  technique  determines  the  three  dimensional  electron 
density  distribution  of  the  ionosphere  over  a  wide  area  at  any  time  from  the  time 
series  of  observed  sounder  data.  The  technique  uses  the  following  approach: 

1 .  Characterize  electron  density  profiles  by  parameter  sets, 

{fs'  ^m/  Ao,  Ai,  A2,  A3,  A4). 

2.  Take  DFT  for  each  parameter  at  each  station  to  determine  the  wave  frequencies 
that  contribute  to  the  observed  electron  density  distribution. 

3.  Determine  the  best  fit  plane  wave  for  each  frequency  from  the  amplitudes  and 
phases  at  the  five  stations. 

4.  Construct  parameter  (e.g.  f^)  contour  map  by  linear  superposition  of  plane 
waves. 

5.  Calculate  electron  density  at  area  point  X,Y  using  parameters  fs,  fm,  Zm/  Aq,  A], 
A2,  A3,  A4,  specified  for  this  point  by  the  constructed  parameter  maps. 

Preliminary  testing  over  the  Northeast  American  region  shows  the  promise  of  this 
technique  in  support  of  OTH  radar,  HF  communications  and  direction  finding 
applications.  The  potential  features  of  the  technique  include  the  capability  of 
providing  maps  in  real  time  for  nowcasting  and  to  make  near-term  forecasts. 

Since  the  technique  is  completely  data  driven,  difficulties  arise  when  data  gaps  in 
the  order  of  hours  or  more  occur  at  any  monitoring  station.  This  problem  requires 
further  studies. 
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